Entropy and correlations in a fluid of hard 
spherocylinders: The onset of nematic and 

smectic order 

D. Costal F. Micali^ F. Saija*, and P. V. Giaquinta'''* 

I Istituto Nazionale per la Fisica della Materia (INFM) and 
Universitd degli Studi di Messina, Dipartimento di Fisica 
Contrada Papardo, CP. 50 - 98166 Messina, Italy 
X CNR - Istituto per i Processi Chimico-Fisici, sez. Messina 
Via La Farina 237 - 98123 Messina, Italy 



Abstract 

Hard spherocylinders (cylinders of length L and diameter D capped 
at both ends with two hemispheres) provide a suitable model for inves- 
tigating entropy-driven, mesophase formations in real colloidal fluids 
that are composed of rigid rodlike molecules. We performed exten- 
sive Monte Carlo simulations of this model fluid for elongations in 
the range 3 < L/D < 5 and up to L/D = 20, in order to investi- 
gate the relative importance of translational and orientational corre- 
lations allowing for the emergence of nematic or smectic order in the 
framework of the so-called residual multi-particle entropy (RMPE). 
The vanishing of this quantity, which includes the re-summed contri- 
butions of all spatial correlations involving more than two particles, 
signals the structural changes which take place, at increasing densi- 
ties, in the isotropic fluid. We found that the ordering thresholds 
detected through the zero-RMPE condition systematically correlate 
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with the corresponding phase-transition points, whatever the nature 
of the higher-density phase coexisting with the isotropic fluid. 

PACS 64.60.-i, 64.70.Md, 61.20.Ja, 65.40.Gr 

1 Introduction 

A systematic approach to the study of the relation between configurational 
entropy, spatial correlations and local ordering in a fluid of elongated particles 
was recently outlined in [1], where the Authors used the theoretical frame- 
work that had been originally introduced by Giaquinta and Giunta in order 
to identify the "hidden" signature of freezing in the structural properties of 
a homogeneous and isotropic fluid of hard spheres [2]. Their approach was 
based on the well known multi-particle correlation expansion of the statisti- 
cal entropy, flrst established for a flnite system by H. S. Green [3] and later 
generahzed to an inflnite open system by R. E. Nettleton and M. S. Green [4]: 

oo 
n=2 

where Sex is the excess entropy per particle in units of the Boltzmann con- 
stant and the "n-body entropies" (s„) are obtained from the integrated con- 
tribution of spatial correlations between n-tuples of particles. Note that the 
one-body term coincides with the entropy of the corresponding ideal gas. The 
pair entropy of a homogeneous and isotropic fluid of nonspherical molecules 
can be written as [5] : 

S2{p) = J {g{r,u;')ln[g{r,u;')]-g{r,u;') + l}drd^^ (2) 

where g{r, o;^) is the pair distribution function (PDF) which depends on the 
relative separation r between a pair of molecules and on the set of Euler angles 
cu^ = [uji, UJ2] that specify the absolute orientations of the two particles in the 
laboratory reference frame. The quantity Q represents the integral over the 
Euler angles of one molecule {Q = 47r for linear molecules), while p is the 
particle number density. 
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The approach originally set up by Giaquinta and Giunta focused on the 
properties of the so-called residual multi-particle entropy (RMPE): 

As = Sex - S2 . (3) 

The quantity As includes, by definition, the re-summed contributions of all 
correlations involving more than two particles [2]. This quantity, despite 
its minor quantitative relevance in the overall entropy balance, is far from 
being irrelevant in that it conceals significant indications on the statistical 
thermodynamics of the fiuid that are intimately related to the role played 
by high-order density correlations in the system. The noticeable feature ex- 
hibited by the RMPE, as compared to the pair entropy (a negative-definite 
quantity), is its non-monotonic behavior. In fact. As (p) starts decreasing 
for increasing densities until the function attains a minimum beyond which 
it undergoes a sharp increase up to positive values. This behavior may come 
out, at first glance, somewhat unexpected in that one would probably sur- 
mise that, in analogy with the monotonic behavior of the two-body term S2, 
the growth of higher-order density correlations that is induced by a further 
compression of the system should throughout contribute to a contraction of 
the configurational states that are accessible to the system for smaller and 
smaller volumes, as is actually the case from low to intermediate densities. 
Instead, one observes a "slowing down" , for increasing pressures, in the reduc- 
tion rate of the average number of states with respect to the non-interacting 
system in an equivalent thermodynamic condition: in fact, the positive cu- 
mulative contribution of multiparticle correlations to the entropic balance 
moderately damps the systematic drop that is caused by the dominant pair 
term. This effect is suggestive of the emergence of a new structural condition 
in the system. More specifically, it indicates that multiparticle correlations 
have started to build up some type of ordering, on a local scale, as a result 
of compeUing entropic constraints. The associated signature is the vanishing 
of the RMPE: a posteriori, it is found to herald in a very sensitive way the 
occurrence - in a proximate range of densities and/or temperatures - of a 
phase transition of the fully disordered fiuid into a more structured phase. 
By now, the reliability of such a one-phase ordering criterion has been exten- 
sively documented for a variety of model fluids [6, 7], also in two dimensions 
[8]. 

With specific reference to a fiuid of hard spherocylinders, i.e., cylinders 
of length L and diameter D capped at both ends with two hemispheres. 
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preliminary evidence was reported in [1] on the RMPE indication concerning 
the isotropic-nematic (I-N) transition. In the case of simple fluids composed 
of spherical particles a number of one-phase criteria have been introduced 
in the past with the aim of locating the transition from the liquid to the 
solid phase without resorting to the knowledge of the free energy of the 
two coexisting phases [9]. However, no such simple criteria are available for 
nematic fluids. In this respect, the zero-RMPE criterion stands as, perhaps, 
a unique exception in that it can be extended in a straightforward way to 
fluids composed of elongated particles to deal with the transition from the 
homogeneous and isotropic fluid to a liquid-crystalline phase. 

Since the seminal work of Onsager [10], the model of hard spherocylin- 
ders has been diffusely investigated, notwithstanding its simplicity, in or- 
der to study excluded- volume effects in real liquid crystals. In fact, hard 
spherocylinders exhibit - depending on the value assumed by the aspect ra- 
tio L/D - a, surprising variety of mesophases [11, 12, 13, 14, 15]. In [1] 
constant-pressure Monte Carlo simulations were carried out for L/D = b 
only, a case where the disordered fluid undergoes, under compression, a flrst- 
order transition to a partially ordered nematic phase. In this paper we have 
extended the above calculations to lower values of the aspect ratio in the 
range 3 < L/D < 5, thus encompassing a region where the homogeneous 
and isotropic fluid undergoes a transition to a fully ordered solid (S) phase 
as well as to the smectic-A (SmA) phase, a layered orientationally ordered 
phase characterized by a one-dimensional density modulation along the direc- 
tion n of the molecular alignment [15] . We have carried out simulations also 
in the limit of very elongated particles, speciflcally for L/D = 20, in order 
to analyze the behavior of the RMPE in a regime where the I-N transition 
point can be determined, asymptotically, in an exact way [10, 16, 17]. This 
is also the range of aspect ratios relevant for the modelling of monodisperse 
aqueous suspensions of rodlike particles like the tobacco mosaic virus, a plant 
virus which shows an isotropic-nematic transition as well as a smectic phase 
at higher concentrations [20]. 

The paper is organized as follows: The theoretical framework and the 
numerical simulation technique, as applied to the case of cylindrical molecules 
with head-tail symmetry, are described in Sects. 2 and 3, respectively. The 
results are presented in Sect. 4 while Sect. 5 is devoted to concluding remarks. 
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2 Theory 



In order to evaluate the RMPE from Eq. (3), one needs to calculate the 

excess and pair entropies of the fluid. The former quantity can be obtained 
by integrating the equation of state (EoS) along a thermodynamic path: 
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In Eq. (4) P is the pressure, [3 is the inverse temperature in units of the 
Boltzmann constant, and p is the density of a suitably chosen reference state. 
The EoS of the fluid was obtained through a Monte Carlo sampling of the 
system at different pressures, while the quantity Scx(p) was calculated using 
the Widom test-particle insertion method [21, 22] at low enough densities. 

As for the pair entropy, this quantity was obtained, according to Eq. (2), 
after integration of the PDF which, in turn, was decoupled in the form [23]: 

g{r,uj^) = g{r)g{uj\) , (5) 

where g{r) is the radial distribution function of the molecular centers of mass 
that is obtained after integrating the full PDF over the set of Euler angles, 
while g{uj'^\r) represents the conditional distribution function, normalized to 
Jl^, for the relative orientation of a pair of molecules whose centers of mass 
lie at a distance r. Using Eq. (5), the pair entropy was separated into the 
sum of a translational (sg"^) and an orientational part (^2'^), after sorting out 
an additional excluded-volume contribution which arises in Eq. (2) from the 
integration over space regions where gi(r, a;^) = 0: 

S2 = -S2P + 4 + S^'-. (6) 

In Eq. (6), B2 is the second virial coefficient of hard spherocylinders: 

S2 = -7rL'^ + 7rL'2L+-L>L^ (7) 
3 4 

The translational contribution is formally analogous to the pair entropy of a 
simple fluid with no rotational degrees of freedom: 

'2--\p l[g(r)\ng(r)-g(r) + l]dr, (8) 



5 



where the integration is carried out over non-overlapping regions. As for the 
orientational term, one finds: 



s^' = P g{r)S°'{r)dr, 



(9) 



where 




(10) 



Even for uniaxial molecules, the computation of the full orientational correla- 
tion function g{u!^\r) and its subsequent integration is a formidable task. In 
this specific case there are only four independent degrees of freedom which 
can be fixed by providing, e.g., the relative intermolecular distance r and 
three angles [24] . Without loss of generality, we can hold the center of mass 
of, say, particle 1 fixed at the origin of the laboratory coordinate system 
and let its molecular axis ei lie along the z-axis. Then, the position and 
orientation of particle 2 can be unambiguously determined by specifying, for 
instance, two polar coordinates of its center of mass, i.e., the distance r and 
the angle formed by the intermolecular axis with the z-a.xis, and the Euler 
angles 9 and of the 62 axis (the "twist" angle, %, is redundant for uni- 
axial molecules). In order to simplify the whole scheme, we neglected the 
dependence of 5^(0;^ |r) on i? and 0, thus focusing on the information associ- 
ated with the angle representing the relative orientation of the two molecular 
axes, 9 = arccos (ei • 62). Indeed, this is the critical parameter that should 
be monitored in order to assess the onset of either nematic or smectic or- 
der in the fiuid. The contracted orientational distribution function g{9\r) is 
normalized in such a way that: 



We do not expect that such a simplification may heavily bias the estimate 
of the weight associated with orientational correlations in the homogeneous 
and isotropic phase. Indeed, the quantity we are ultimately interested in is 
just an integral measure that can help to assess the relative balance between 
the positional and angular contributions that intervene in the configurational 
entropy. 
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3 Simulation 



We performed Monte Carlo (MC) simulations of hard spherocylinders at 

constant pressure. This choice is rather common for systems composed of 
non-spherical hard-core particles for which it may be difficult to calculate 
the equation of state through the contact value of the PDF, as is usually 
done in a constant-volume simulation. Furthermore, density fluctuations 
are more naturally accommodated in a constant-pressure simulation [22]. 
The number of particles, that were enclosed in a cubic box with periodic 
boundary conditions, typically ranged between 500 {L/D = 3.2,5) and 800 
(L/ D = 3), with an intermediate value of 672 particles for L/D = 4. A larger 
sample composed of 1500 particles was used for more elongated particles with 
L/D = 20. 

A perfectly ordered, low-density lattice arrangement was the starting con- 
figuration. Then, the system was allowed to equilibrate, initially along a 
const ant- volume MC path. A sequence of states was then generated at con- 
stant pressure, starting from an equilibrated isotropic configuration. We 
carried out rather extended simulations, as compared with those currently 
reported in the literature, especially in proximity of the transition threshold. 
In fact, after an equilibration period of typically 10^ to 10^ MC cycles, simu- 
lation data were obtained by generating chains consisting of 5 x 10^ to 5 x 10^ 
MC cycles, depending on the pressure. A MC cycle consisted of an attempt of 
changing, sequentially, both the position and the orientation of each molecule, 
followed by an attempt of changing the volume of the sample. During such 
a procedure, each molecule was first displaced from its original position, by 
moving at random its center of mass within a cube. The molecule was then 
re-oriented within a cone centered about the molecular axis. The overall ac- 
ceptance probabihty for the combined move was adjusted to about 50%. The 
isotropic change of the simulation box was obtained by varying at random 
the box side L. Finally, the maximum displacement allowed was set so as to 
achieve an average acceptance probability for volume changes of ~ 50%. For 
hard-core molecules, the usual Metropolis acceptance criterion reduces to a 
check of particle overlaps after each Monte Carlo move. In this regard, a sim- 
ple and efficient algorithm can be used for spherocylinders [25, 26]. During 
the cumulation runs, we monitored the average density and built histograms 
for the translational and orientational distribution functions, respectively; 
equilibrium averages and standard deviations were computed by dividing 
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such runs into 5 to 10 independent blocks. In the following, we shall refer 
to the volume packing fraction r] = pVhsc, where Vhsc = (7r/4)Z}^L + {7i/6)D^ 
is the volume occupied by one spherocylinder, and to the reduced pressure 
P* = pPVhsc- We sampled the radial distribution function g{r) at intervals 
of O.ID, while the orientational distribution g{9\r) was calculated with a grid 
spacing of A^^ = 5°. 

We investigated the phase diagram of the model up to the threshold of 
either the nematic or the smectic phase. The onset of hquid-crystalline order 
was monitored through the nematic order parameter: 



where ■0i is the angle formed by the molecular axis of the i-th molecule 
(the sum running over all particles in the system). The director n defines 
the alignment direction attained, on average, by the molecules. Since this 
direction is not known a priori, the order parameter is usually determined by 
diagonalizing the traceless, symmetric, second-rank tensor Q that is defined 
in terms of the molecular axes as: 



where I is the unit tensor. The nematic order parameter is then obtained as 
the largest eigenvalue, and the director as the corresponding eigenvector [25]. 
The value attained by S is close to zero in the isotropic phase and tends to one 
in a highly ordered phase. Across the ordering transition, both the density 
and the order parameter show a finite jump. However, such discontinuities, 
that are associated with a first-order phase transition, are in certain cases so 
small that a closer investigation of the correlation functions or a direct obser- 
vation of the arrangement of the molecules in the sample may be necessary 
on the thermodynamic edge of both the nematic and smectic phase. 

As already anticipated, we used Widom's insertion method [21, 22] to ob- 
tain independent numerical estimates of the excess entropy at low densities, 
=^cx(p). to be used in Eq. (4). The Widom technique provides a simple scheme 
for calculating the excess chemical potential, ficx, of not too dense a liquid via 
the calculation of the average Boltzmann factor associated with the (never 
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accepted) random insertion of an additional particle in an A^-particle sys- 
tem. The excess entropy can then be obtained through the thermodynamic 
equation: 

Sex = H 1 (14) 

P 

Depending on the density and the elongation of the spherocylindcrs, 50 to 
200 trial insertions per MC step were typically needed in order to obtain a 
reliable estimate of /igx- The use of this technique at small but finite densities 
avoids the otherwise necessary extrapolation to the ideal gas limit, as a refer- 
ence state for the integration contemplated in Eq. (4). This latter procedure 
was adopted in [1], where some spurious fluctuations of the orientational dis- 
tribution function as well as some difficulties associated with the simulation 
in the highly dilute regime were registered. 

As sketched in Fig. 1, the values obtained for Sex using Widom's insertion 
method arc quite stable in all the cases we have examined: the thermody- 
namic potential relaxes almost immediately, with small fluctuations about 
the average value. No significant drifts were observed during simulations as 
long as 10^ MC steps. The simulation parameters and the results for the 
relevant properties of the fiuid at low densities were collected for various 
elongations in Table 1. 



4 Results and discussion 

4.1 Phase behavior for 3 < L/ D < 5 

In this subsection we shall focus on the phase behavior of the model in a range 
of parameters where, for increasing values of the aspect ratio, the disordered 
fluid undergoes a transition to either a solid, a smectic or a nematic phase. 
We investigated the model for shape anisotropics L/D = 3,3.2,4, and 5. 
For L/D — 3, the fluid freezes into a crystalline phase. An I-SmA phase 
transition is observed over a range of larger aspect ratios bounded by an 
isotropic-smectic-solid (I-SmA-S) triple point {L/D — 3.1) and by an I-N- 
SmA triple point [L/D = 3.7), respectively. For even larger anisotropies, the 
disordered fluid forms, under compression, a nematic phase [15]. 

The rehability of the present results for the EoS and for the nematic 
order parameter can be appraised through Fig. 2 where the above quantities 
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were plotted as a function of the packing fraction for L/D = 3.2 and 5. 
The comparison with the simulation data reported by McGrother and co- 
workers [14] shows that a system with 500 particles is fairly adequate for 
evaluating the basic thermodynamic properties of the fluid. The agreement 
turned out to be equally good along the disordered branch also for LjD — 2> 
and 4, up to the transition point. 

We now turn to a discussion of the properties of the model for L/Z^ = 5, a 
case where the isotropic fluid undergoes a rather weak first-order transition 
to a thcrmodynamically stable nematic phase at about 45% of the close- 
packing density, a value corresponding to an absolute packing fraction of 
0.4 [14, 15]. No specific feature indicating the transition edge is apparent 
either in the equation of state or in the order parameter (see Fig. 2). A 
rather modest degree of positional order shows up in the spatial profile of 
the radial distribution function that was plotted in Fig. 3 for increasing 
pressures, all the way up to the transition threshold. Actually, for moderately 
small packing fractions, the local density around one generic sphero cylinder 
is significantly reduced, with respect to the homogeneous value taken at large 
distances, for interparticle separations lower than ~ 3| molecular diameters, 
a distance which corresponds to the minimum contact separation between 
two orthogonal spherocylinders. 

Retrospectively, this local "depletion" effect is more easily understood if 
one looks at the orientational distribution function (see Fig. 4): the almost 
fiat profile of 5'(^|r) for r > 3.5D, when plotted as a function of ^, implies that 
orthogonal geometries in the relative spatial arrangement of two spherocylin- 
ders arc equally probable as parallel configurations are. The steric hindrance 
produced, on average, by such orthogonal (or almost orthogonal) particles 
sensitively reduces the probability of finding other molecules at shorter dis- 
tances from the central one. Upon further compressing the fluid, the value 
of g{0\T) for aligned-particles geometries {§ — 0, tt) steadily increases while 
different geometries become, at the same time, rarer and rarer. This is ap- 
parent from the progressive and systematic reduction of g(Q\r\ for short 
and intermediate molecular separations, over a wide range of angles centered 
about 7r/2. Correspondingly, the radial distribution function becomes more 
and more structured at short distances as is distinctly shown in Fig. 3 by the 
building up of a first coordination shell formed by particles lying almost at 
close contact with the central one. 

The excess entropy of the fiuid was resolved in Fig. 5 in terms of the 
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contributions associated with pair and more-than-two-particle correlations, 
respectively (see Eq. (3)). The behavior of the RMPE is similar to that 
exhibited by the same quantity in a variety of other model systems, composed 
of spherical particles, across an ordering transition such as the freezing of a 
fluid [6, 8] or the phase separation of a binary mixture [7]. In particular, 
this quantity vanishes for a value of the packing fraction {rjo) which falls 
precisely on top of the currently estimated I-N phase transition point (see 
Table 2 where we summarized the predictions of the entropy-based criterion 
for the presently investigated values of the aspect ratio). We note that the 
current RMPE estimate sensitively improves on that formerly given in [1] 
where the values reported for the pair entropy turned out to be affected 
by a systematic computational error. However, no significant morphological 
change resulted, a posteriori, in both the pair entropy and in the RMPE. We 
also recall that the Authors had shown in [1] that the pair entropy changes 
with the size of the sample, for a given density, in a tiny but systematic 
way, the zero-RMPE density shifting towards moderately larger values. The 
RMPE inverts its trend (from decreasing to increasing) in the density region 
(77 ~ 0.36) corresponding to a more extended structuring of the fluid at short 
distances (see Fig. 3). To be specific, we refer to the growing of a local 
nucleus of (mostly) aligned particles that is signalled by the emergence of a 
secondary maximum in the radial distribution function at a relative distance 
twice as large as that corresponding to the first coordination shell. 

Figure 5 also shows the pair entropy of the fluid resolved, according to 
Eq. (6), into the sum of a (spatially independent) excluded- volume term 
(— i?2p) plus two more contributions arising from translational and orienta- 
tional correlations, respectively. As expected, it is s™' that drives the rather 
sharp increase of As{r]) just beyond the minimum. Vice versa, the contri- 
bution associated with positional correlations, s*2j ^^^Y small all over the 
fluid range. The function which, upon integration (see Eq. (9)), yields the 
quantity is shown, for increasing densities, in the right part of Fig. 3. 
The value of the corresponding integral is negative and keeps close to zero at 
low and intermediate densities (see Fig. 5). For packing fractions larger than 
~ 0.3, starts decreasing more and more rapidly and eventually overcomes 
(in absolute value) the excluded-volume term. As seen from Fig. 5, the fln- 
gerprint of the ordering of the fluid is manifest in the persistence of angular 
correlations up to intermediate and long distances. 

Considerations similar to those developed above for L/D = 5 apply also 
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to the L/D = 3.2 case, even if it is a smectic phase that is now formed by 
the isotropic fluid. However, the basic underlying mechanism which drives 
the first-order I-SmA transition is analogous to that exploited in the nematic 
case. Obviously, the denser isotropic fluid looks more structured than for 
L/D = (see Fig. 6). Furthermore, at the transition point angular correla- 
tions show a very slow spatial decay. We also note that in this case the lower 
molecular asymmetry modifies the scale of interparticle distances in the av- 
erage density profile. Thus, the effect associated with the structuring of the 
second coordination shell, which we already commented upon for L/D = 5, is 
less resolved since the corresponding distance falls very close to the position 
of the lone weak maximum, located at a distance r/D = 1 + ^{L/ D) = 2.6, 
that is present in the radial distribution function at low densities. A shoul- 
der, rather than a peak, emerges at r/D ~ 2.4, again in a density range 
corresponding to the position of the minimum [r] ~ 0.47) attained by the 
RMPE (see Fig. 7). 

In general, the ordering thresholds detected through the vanishing of the 
RMPE correlate very well with the corresponding transition point (see Ta- 
ble 2), whatever the nature of the higher-density phase coexisting with the 
isotropic fluid. The evidence discussed so far justifies the claim that the 
RMPE of sphero cylinders monitors the leading microscopic process which 
drives the phase transition in all such cases, i.e., the aligning of sphero- 
cylinders along a common direction. As such, the vanishing of this quantity 
presumably portends the density threshold beyond which the nematic or- 
dering may possibly emerge in the system, regardless of other concurrent 
structural phenomena like those which may ultimately lead to the formation 
of a more ordered macroscopic state, be it smectic or crystalline, instead of 
a pure nematic phase. An indication consistent with this thesis that is, per- 
haps, not accidental (yet to be considered with due caution because of the 
numerical uncertainty arising from the computational error which affects, at 
least, the third significant figure) is offered by the apparently greater accu- 
racy of the predictions given by the RMPE-based criterion for L/D = 4 and 
5 (see Table 2), viz., in those cases where the stable coexisting phase is a 
truly nematic fluid. In the other two cases investigated, i.e., when either a 
solid or a smectic phase is formed, the criterion just overshoots the target, 
locating the internal threshold for the ordering of the fluid at a shghtly higher 
density. 
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4.2 The Onsager regime 

As the aspect ratio increases, the I-N transition point moves to lower and 
lower densities. In the Onsager limit of infinitely long rods {L/D — > oo), 
a fluid of hard sphero cylinders undergoes a flrst-order phase transition to a 
nematic phase for B2P — 3.29 [10, 16, 17]. However, the phase behavior char- 
acteristic of the Onsager regime qualitatively sets in for shape anisotropies 
even less than L/D = 20, even if for such values of the aspect ratio the basic 
assumption made by Onsager that all virial coefficients of order higher than 
two can be neglected is not yet satisfied [15]. In fact, the contribution of the 
third virial term to the free energy, evaluated for L/D = 20 at the transition 
density, turns out to be still comparable to the contribution of the second 
virial term [18, 19]. 

We carried out simulations for L/D = 20. In passing, we recall that 
this value mimics the typical size of the tobacco mosaic virus, whose parti- 
cles have a practically rigid rod shape of 15 — 18 nm diameter and 300 nm 
length [20]. For such large values of the aspect ratio, the simulation of a 
fiuid of hard spherocylinders dictates more severe constraints on the size of 
the simulation box which should be large enough to accommodate at least 
two rod lengths, so as to avoid multiple overlaps between a molecule and 
the periodic images of another particle [15]. The number of particles cho- 
sen (A^ = 1500) is sufficient to ensure that the condition mentioned above 
is fulfilled in the (relatively low) density range where the disordered phase 
is thermodynamically stable. The reduced pressure and the nematic order 
parameter are shown in Fig. 8. The isotropic liquid phase is mechanically 
stable up to a reduced pressure P* — 0.92: in this regime the system shows, 
for simulations as long as 2.5 x 10^ MC cycles, no significant fluctuations 
about the average values of the packing fraction ((77) = 0.145) and of the 
nematic order parameter {{S) = 0.058), respectively. Upon further com- 
pressing the system up to P* = 0.93, we observed, after about 1.5 x 10^ MC 
cycles, a spontaneous and rapid transformation into a nematic phase which, 
at equilibrium, settled down at a packing fraction (r)) — 0.182 and was char- 
acterized by an order parameter (S) = 0.862. The resulting "boundaries" are 
slightly shifted toward higher densities if compared with the direct estimates 
reported in [15], i.e., r/iso = 0.139 or 0.137, and ?7nem = 0.171 or 0.172, accord- 
ing to whether the Gibbs-Duhem integration or a Gibbs-ensemble technique 
was used. Exploratory runs performed by decompressing the system from 
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the ordered phase distinctly show an hysteresis, the fluid remaining almost 
nematic for pressures as low as P* ~ 0.75. We calculated the pressure also 
using the "decoupling approximation" between the translational and orien- 
tational degrees of freedom originally introduced by Parsons [27]. In this 
approximation a system of rods is mapped onto a system of spheres interact- 
ing with an effective orientation-dependent potential. As already observed by 
McGrothcr and co-workers for the shorter elongation regime 3.2 < L/ D < 5 
that was analyzed in [14] , the resulting EoS turns out to be fairly accurate in 
the isotropic phase up to the transition point, and definitely better than the 
prediction based on the scaled-particlc-theory approach exploited by Boub- 
lik [28]. However, it should be noted that the latter approach is based on 
a one-parameter approximation for the third and fourth virial coefficients 
which cannot be expected to yield accurate results for large elongations [29] . 
A better agreement with the simulation data is obtained with a closely re- 
lated description that was proposed by Nezbeda [29, 30]: this approximation 
is based on an empirical fit of the virial coefficients, especially tailored on 
the available results for hard sphero cylinders. All the above approximations 
for the EoS were reported, for a comparison, in Fig. 8: note that the simula- 
tion data are closely bracketed by the Parsons and Nezbeda approximations, 
respectively. 

Radial correlations are shown in Fig. 9 together with the function whose 
space integral yields sf. All over the liquid density range, the radial distri- 
bution function is practically structureless apart from the correlation hole at 
short distances (see Fig. 9). As a result, the translational contribution to 
the pair entropy keeps close to zero all over the explored density range (see 
Fig. 10). The same is also true for the orientational contribution for packing 
fractions less than about 0.15. It thus follows that the pair entropy of the 
isotropic fluid is substantially equal to the excluded-volume term —B2P, as 
is actually assumed in the original Onsager theory. At variance with the 
positional term, the angular contribution is responsible, together with the 
excluded-volume term, of the vanishing of the RMPE. Just beyond the tran- 
sition point, 52^^ blows up and deflnitely overcomes the excluded-volume term 
(see Fig. 10). The resulting RMPE shows the typical behavior already ob- 
served for smaller aspect ratios. Also in this case, the zero-RMPE threshold 
is in very good agreement with the I-N transition point estimates given in 
[15] (see Table 2). We also reported in Fig. 10 the Parsons and Nezbeda 
analytical estimates for the excess entropy. As already pointed out for the 
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EoS, both approximations turn out to predict fairly well, almost with the 
same degree of accuracy, also the excess entropy of the fluid. 

As seen from Figs. 5, 7 and 10, the overall importance of the cumulative 
contribution of higher-order density correlations to the excess entropy de- 
creases with increasing shape anisotropies. This is clearly a consequence of 
the progressive confinement of the isotropic phase to lower and lower densi- 
ties. In this regard, we just observe that the value attained by the RMPE at 
its minimum drops from ~ 45% of the excess entropy ior L/D — 3.2 to ~ 18% 
for L/D = 20. This behavior is consistent with the numerical findings on 
the virial expansion cited above and with the asymptotic trend predicted by 
Onsagcr as L/D ^ oo. Note, however, that at low densities the ratio As/scx 
moderately increases with L/D for a given packing fraction. Actually, longer 
and thinner particles give rise to more extended multiparticle correlations 
which trigger the ordering of the fluid at lower and lower densities. Beyond 
the zero of the RMPE, where the nematic phase is thermodynamically stable, 
the present data monitor de facto a system in a metastable, macroscopically 
isotropic condition whose proper description lies beyond the standard On- 
sager theory. 

5 Concluding remarks 

In this paper we have analyzed the contributions to the "pair entropy" S2 of 
hard spherocylinders of length L and diameter D that are associated with 
excluded- volume, translational and angular correlations, with particular em- 
phasis on the thermodynamic behavior exhibited by the above quantities on 
approaching the phase transition from the low-density fully disordered phase 
to a partially ordered phase at higher densities. An elucidating synthesis of 
the results obtained for S2 through an intensive Monte Carlo sampling of the 
model is offered by Fig. 11 where we plotted the three contributions cited 
above, each of them being calculated in the state corresponding to the van- 
ishing of the residual multi-particle entropy (RMPE) for an assigned value of 
L/D. The RMPE is a quantity that gauges the cumulative weight of spatial 
correlations involving more than two particles in the configurational entropy 
of the fluid. Figure 11 neatly indicates the relative importance of the three 
contributions in the overall entropic balance: for very small asymmetries, 
the translational term, S2, prevails (in absolute value) while being compa- 
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rable with the excluded- volume contribution, —B2P (note that the data for 
L/D = do obviously refer to the hard-sphere case [2]). For increasing 
values of the aspect ratio the angular term, s^"^, acquires more and more im- 
portance: in fact, an interpolation of the available data shows that the three 
contributions reach the same magnitude for L/D c:^ 0.4. Actually, it is in 
this range of values (more precisely, for L/D < 0.35 ± 0.05 [15]) that the 
"rotator phase" (i.e., the plastic crystal) disappears in the phase diagram of 
hard spherocylinders, and is replaced at high densities by an orientationally 
ordered, crystal phase. As seen from Fig. 11, in the range 0.4 < L/D < 5 the 
ordering of the fluid is mainly controlled by angular correlations. In fact, the 
translational contribution rapidly tends to smaller and smaller values and 
practically vanishes for L/D > 5. Also ^2^^ asymptotically tends to zero in 
this range of increasing asymmetries but, as seen from Fig. 11, the decay 
of this quantity is much slower than that shown by its translational coun- 
terpart. Correspondingly, the excluded-volume term becomes the dominant 
contribution to the configurational entropy of the fluid. Note that its value 
for L/D = 20 (—3.45) is already very close to the Onsager limit (—3.29). 

As far as the quahty of the "predictions" offered by the zero-RMPE crite- 
rion are concerned, the agreement with the properly defined thermodynamic 
thresholds is impressive for all the asymmetries that have been investigated. 
Indeed, the vanishing of the RMPE monitors, in a very sensitive and reli- 
able way, the ordering of the homogeneous and isotropic fluid into a more 
ordered phase be it nematic, smectic or solid. However, we remark that such 
a close quantitative correspondence should not be necessarily and systemat- 
ically expected a priori. In fact, the zero of As needs not coincide with the 
thermodynamic edge of the disordered phase that is properly located through 
the comparison of the free energies of the two coexisting phases. Actually, the 
RMPE conveys a type of information that is intermediate between the fully 
macroscopic level and the underlying microscopic description of the system. 
The interpretation of such a simply accessible information (and, more specif- 
ically, of its change of sign) as an intrinsic signature of the incipient ordering 
of the fluid appears to be quite flrmly established on the basis of an ample 
and coherent evidence that has, by now, emerged in a variety of continuum 
and lattice model fluids [6, 7, 8]. The present analysis of the "fine structure" 
of the configurational entropy of hard spherocylinders further corroborates 
the sensitivity of the RMPE to local ordering phenomena which may even 
turn out into a rather weak first-order phase transition. We believe that this 
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one-phase entropic criterion can be conveniently generalized to account for 
more than one phase transition, leading the system to increasing levels of 
self-organization. An attempt in this direction has been recently done by 
Cuetos and coworkers whose results seem to indicate the existence of a jump 
in the translational part of the pair entropy as well as in its rate of change 
across the nematic-smectic transition [32]. However, this evidence was not 
considered conclusive because of the limited size of the simulation box. We 
are also currently investigating this point with promising preliminary results. 
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Tables 



Table 1: Thermodynamic properties of the fluid in the low-density states p 
used for the integration of the equation of state in Eq. (4) of the text. Column 
VllI refers to the number of trial insertions per MC cycle; the percentage ratio 
of accepted trial insertions is reported in column IX. Standard deviations in 
the last digit (s) are given in parenthesis. 



L/D P* [f]) (sex) (S2) MC cycles Trials Ratio (%) 



3 


0.30 


0.131(1) 


2.33(1) 


-1.039(13) 


-0.959(2) 


5 X 10^ 


100 


9.7 


3.2 


0.30 


0.130(1) 


2.37(1) 


-1.054(17) 


-0.971(3) 


9 X 10^ 


50 ^ 200 


9.4 


4 


0.30 


0.125(1) 


2.55(1) 


-1.149(6) 


-0.996(1) 


5 X 10^ 


80 200 


7.8 


5 


0.181 


0.089(1) 


1.90(1) 


-0.878(6) 


-0.834(1) 


6 X 10^ 


20 


14.9 




0.379 


0.136(1) 


3.25(1) 


-1.457(24) 


-1.237(4) 


5 X 10^ 


200 


3.8 


20 


0.10 


0.044(1) 


2.39(1) 


-1.137(13) 


-1.066(1) 


5 X 10^ 


50 


9.2 




0.20 


0.066(1) 


3.79(1) 


-1.785(26) 


-1.588(5) 


10 X 10^ 


200 


2.2 



20 



Table 2: The packing fractions corresponding to the vanishing of the residual 
multiparticle entropy are compared, for different elongations, with the ther- 
modynamic thresholds of the phase transitions undergone by the isotropic 
fluid, according to the available simulation data. Note that the value as- 
cribed to Ref. [15] for L/D = 4 was estimated through a hnear interpolation 
of the data for the I-N transition densities reported in Ref. [15] for L/D — 3.8 
and L/ D = 4.2, respectively. A slightly lower estimate for the I-N transition 
density L/D — A, r] — 0.459, is reported in Ref. [31]. 



Phase transition point 



L/D Phase transition 770 Ref- [14] Ref. [15] 



3 

3.2 

4 
5 

20 



I-S 
I-SmA 
1-N 
1-N 
I-N 



0.521 
0.510 
0.463 
0.399 
0.148 



0.490 
0.513 
0.472 
0.407 



0.512 
0.503 
0.466 
0.398 
0.139 
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Figures 



Figure 1: Typical evolution of Sex(p) (see Table 1) during a MC run for 
different elongations. The data were obtained using the Widom test method. 
Circles: L/D = 3.2 at P* = 0.30; diamonds: L/D = 5 at P* = 0.379; 
squares: L/D = 20 at P* = 0.20. 



Figure 2: Equation of state (top) and nematic order parameter (bottom) 
for L/D — 3.2 (left panels) and 5 (right panels). Solid circles: this work; 
triangles: MC simulations by McGrother et al. [14]. The error bars are 
systematically smaller than the size of the markers. 



Figure 3: Radial distribution function (left) and integrand (right) for s'^ 
in Eq. (9) at several densities for L/D = 5. Open triangles: P* = 0.18, 
{n) = 0.089; solid circles: P* = 1.10, (77) = 0.223; squares: P* = 2.53, 
[r]) = 0.310; solid triangles: P* = 4.20, (77) = 0.372; diamonds: P* = 4.95, 
(ri) = 0.397. 
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Figure 4: Orientational distribution functions g{9\r) at P* — 1.10 (left) and 
4.95 (right), plotted, from top to bottom, for several distances {r/D — 1.25, 
1.95, 2.55, 3.65, 6.05). 

Figure 5: Left panel: residual multiparticle entropy (circles) resolved into 
the excess (squares) and pair (solid diamonds) contributions for L/D — 5; 
the values of the excess entropy at lower densities were also computed us- 
ing the Widom test method (pluses). Right panel: pair entropy (solid di- 
amonds) resolved into the translational (triangles), orientational (squares), 
and excluded- volume (circles) contributions, according to Eq. (6). Lines are 
smooth interpolations of the simulation data. 

Figure 6: Radial distribution function (left) and integrand (right) for s'^ in 
Eq. (9) at several densities for L/D — 3.2. Open triangles: P* — 0.5, {rf) = 
0.170; solid circles: P* = 2.0, (77) = 0.305; squares: P* = 5.0, (r?) = 0.411; 
solid triangles: P* = 8.0, {r]) = 0.474; diamonds: P* = 9.8, (r?) = 0.507; 
pluses (right panel): P* = 9.9, {rj) = 0.511. 

Figure 7: Left panel: residual multiparticle entropy (circles) resolved into the 
excess (squares) and pair (solid diamonds) contributions for L/D = 3.2; the 
excess entropy at the lowest packing fraction was computed using the Widom 
test method. Right panel: pair entropy (solid diamonds) resolved into the 
translational (triangles), orientational (squares), and excluded- volume (cir- 
cles) contributions, according to Eq. (6). Lines are smooth interpolations of 
the simulation data. 

Figure 8: Left panel: equation of state for L/D = 20 (solid circles); the 
Parsons (dotted line, [27]), Boublik (solid line, [28]), and Nezbeda (dashed 
line, [30]) approximations for the equation of state are also shown. Right 
panel: nematic order parameter plotted as a function of (?]). Triangles refer 
in both panels to states obtained upon decompressing the nematic fluid. The 
error bars are systematically smaller than the size of the markers. 



23 



Figure 9: Radial distribution function (left) and integrand (right) for in 
Eq. (9) at several densities for L/D = 20. Circles: P* = 0.10, {r]) = 0.044; 
squares: P* = 0.50, {r]) = 0.107; triangles: P* = 0.85, {r]) = 0.140; solid 
diamonds: P* = 0.93, [r]) = 0.182; pluses (right panel): P* = 0.70, {r]) = 
0.127. 



Figure 10: Left panel: residual multiparticle entropy (circles) resolved into 
the excess (squares) and pair (sohd diamonds) contributions for L/D — 20; 
the excess entropy at lower densities was also computed using the Widom 
test method (pluses). The Parsons (dotted line, [27]) and Nezbeda (dashed 
line, [30]) approximations for the excess entropy are also shown. Right panel: 
pair entropy (solid diamonds) resolved into the translational (triangles) , ori- 
entational (squares), and excluded- volume (circles) contributions, according 
to Eq. (6). Solid hues are smooth interpolations of the simulation data. 



Figure 11: Translational (circles), orientational (squares), and excluded- 
volume (solid diamonds) contributions to the pair entropy of the fluid, cal- 
culated in the states of vanishing residual multiparticle entropy, for different 
shape anisotropics. Data for L = (hard spheres) were taken from Ref. [2]; 
data for L = 1 were obtained from the equation of state reported in Ref. [13]. 
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